3-D transformations. ---------------------- "Several large, artificial constructions are approaching us," ZORAC announced after a short pause. "The designs are not familiar, but they are obviously the products of intelligence. Implications: we have been intercepted deliberately by a means unknown, for a purpose unknown, and transferred to a place unknown by a form of intelligence unknown. Apart from the unknowns, everything is obvious." -- James P. Hogan, "Giants Star" 3D transformations would assume the central position in the engine. After all, this is something which would allow us to see objects from different sides, manage the three dimensional space etc. The simplest but by no means least important transformations are translation and scaling. The rotation transformation would be little bit more complicated to understand and implement but fortunately not much. Talking of all three transformations we can either consider them from the point of an object (collection of points in 3D) being transformed, or a space itself being transformed. In some cases it is more productive to think one way in some another way. Translation and scaling. ------------------------ Translation. This transformation moves points of an object to a new specified position pic 2.1, On the other hand we may think of it as a system of references being moved in the opposite direction pic 2.2: Y ^ ^ | * A' | ^ * B' | | | ^ Y'^ |A*--------+ | y component | |A* | B*--------+ | B* | x component | | ---------+-------------------> X - - - - - - | - - - - - - -> X 0| | 0| | -----------+------------->X' | 0'| | | | | | | | pic 2.1 pic 2.2 This way if it is a collection of points we are moving, Distances between them would not change. Any move across 2D space can separate into 2 components - move along X and move along Y. It would take 3 components in 3D space - along X, Y and Z. Important place we may need translation is to , for example , move (0,0) of the screen from upper left corner of the screen to it's physical centre. So the translation function might look like: void T_translation(register int *from,register int *into,int number, int addx,int addy,int addz ) { register int i; for(i=0;i X 0| | | pic 2.3 The obvious usage of this transformation: trying to fit 2000x2000 object into 200x200 window - evidently every point of this object would have to have coordinates scaled to fit in the window range. In this case scaling is just multiplication by 1/10. The object is contracted that way. On the other hand we can think that it was actually the window space expanding to enclose the object (we would have to abstract far enough from computer hardware today on the market to do that, however). In general case we can introduce separate factors for every axis, if this factor is greater then 0 and less then 1 we would contract the object if it would be greater then 1 the object would be expanded. The code is once again trivial: void T_scaling(register int *from,register int *to,int number, float sx,float sy,float sz ) { register int i; for(i=0;iX | x | pic 2.5 Now, let's consider the situation when we rotate the system of references counterclockwise by the angle alp. What would be coordinates of point A in the turned system of references Y'X' if it's coordinates in the original system of references were (x,y)? ^ Y X' | / Y' Y/*\--*A / \ / | \| / \ / | |\ / Yy'\ | | /Yx' \ | | \| /\ Xx' -------------+---+-------------------> X /| \/ X / Xy' \ / | \ / | \ / | \ | pic 2.6 Using the sin/cos formulas from above we can find projections of x and y (that is of the projections of the point onto original axis's) onto new axis's X'and Y'. Let's sum projection of x onto X' and projection of y onto same axis's and do the same thing with projections of x and y onto Y'. Those sums would be just what we are looking for, coordinates of the point A in the new system of references X'Y' i.e. X' = Yx'+Xx' = y*sin(alp) + x*cos(alp) Y' = Yy'+Xy' = y*cos(alp) + (-x*sin(alp)) Note the sign of Xy', just from looking on the pic2.6 one can see that in this case x is being projected onto the negative side of Y'. That's why the negative sign. So those are the transformation formulas for counterclockwise rotating system of reference by an angle alp: x'= y*sin(alp) + x*cos(alp) y'= y*cos(alp) - x*sin(alp) On the other hand we can consider it as point A itself rotating around zero clockwise. 3-D rotation. ------------- This subject is complicated only because it is a bit hard to visualize what is happening. The idea on the other hand is simple: applying three 2-D rotations one after another so that each next works with the results obtained on the previous stage. This way every time we are applying new transformation we are working not with the original object but the object already transformed by all previous rotations. Each of 2D rotations this way is bound to some axis around which two other axes would rotate. There are few important considerations influencing the way we would find the formulas: 1) What kind of system of references we have. What is the positive direction of the rotations we would be applying; and 2) In what order we want to apply 2D rotations. System of references: --------------------- Well, we do have freedom of choosing directions of axes, for one thing, we are working in orthogonal system of references where each axes is orthogonal to the plane composed of remaining two. As to where the positive direction of every axis point we have freedom to decide. It is customary for example to have positive side of Y directed upside. We don't actually have to follow customs if we really don't want to, and remembering that the colourmap, we would be doing all drawings into, have Y directed down we might want to choose corresponding direction here. However having it pointing up is more natural for us humans, so we might save time on debugging a while afterwards. Let's point it up. As to the Z let's direct it away from us and X to the right. Y^ Z | / | / alp /|<---+ | |/ | -------|-+-----------> X bet | | __ V/| /| gam /----+ / | / | | pic 2.7 As to the rotations, let's call the angle to turn XY plane around Z axis as gam, ZY around X as bet and ZX around Y as gam, pic 2.7 Transformation order. --------------------- The problem with transformation order is that the point consequently turned by gam-bet-alp won't be necessarily at the same place with where it might go if turned alp-bet-gam. The reason is our original assumption: each next transformation works on already transformed point by previous 2D rotations. That is, we are rotating coordinates around moving axes. On the picture 2.7 we can see that rotation bet is performed around axis x but the next rotation alp is performed around axis z which after application of previous transformation would move. + +----+ z / /| / / alp / + V bet \ ------------+-------|---- x ---------+------------ x / | \ +----+ z | / | | V alp pic 2.8 The implication is: we have to know with respect to what each next rotation is performed, that is what is the order of applications of 2D rotations. Let's think how we are coordinating the surrounding world? First we think of the direction. Then we can tilt our head up and down, and finally from there we can shake it from left to right. Now, when we are tilting head we have already chosen the direction. If we first would tilt head then the directional rotation to be performed after would not be parallel to the ground but rather perpendicular to the imaginary line being continuation of our tilted neck. ( No responsibility hereby assumed for all direct or consequential damage or injury resulted from doing that ) So it all comes to directional rotation first - gam in our case. pitch second - bet; roll last - alp. Let's do just that using obtained formula for 2D rotation: ^Z ^Y' ^Y" | | | |<---+ |<---+ |<---+ | gam| | bet| | alp| | | | | | | -+----+--->X -+----+--->Z' -+----+--->X" | | | x'=z*sin(gam)+x*cos(gam) y'=y z'=z*cos(gam)-x*sin(gam) x"=x y"=y*cos(bet)-z*sin(bet) z"=y*sin(bet)+z*cos(bet) x"'=y*sin(alp)+x*cos(alp) y"'=y*cos(alp)-x*sin(alp) z"'=z pic 2.9 That's basically it, using those formulas we can apply 3D rotation to coordinates x,y,z and get x"',y"',z"'. We can see here 12 multiplication's, The question is can we reduce this number? Let's try getting rid of temporary variables x',y',z',x",y",z". We can do it by just substituting each occurrence of them by their real meaning: First let's try obtaining x",y" and z" expressed via x,y,z: y"=y'*cos(bet) - z'*sin(bet) y"=y*cos(bet) - (z*cos(gam)-x*sin(gam)) *sin(bet) y"=y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - x"=x' x"= z*sin(gam) + x*cos(gam) - - - - - - - - - - - - - - - z"=y'*sin(bet) + z'*cos(bet) z"=y*sin(bet) + (z*cos(gam)-x*sin(gam)) *cos(bet) z"=y*sin(bet) + z*cos(gam)*cos(bet) - x*sin(gam)*cos(bet) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Now using that we can express x"',y"' and z"' via x,y,z: y"'=y"*cos(alp) - x"*sin(alp) y"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*cos(alp) - (z*sin(gam) + x*cos(gam)) *sin(alp)) y"'=y*cos(bet)*cos(alp)-z*cos(gam)*sin(bet)*cos(alp)+x*sin(gam)*sin(bat)*cos(alp) -z*sin(gam)*sin(alp) - x*cos(gam)*sin(alp) y"'=x* [ sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)] + y* [ cos(bet)*cos(alp) ] + z* [-cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)] ----------------------------------------------------------- x"'= y"*sin(alp) + x"*cos(alp) x"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*sin(alp) + (z*sin(gam) + x*cos(gam)) *cos(alp) x"'=y*cos(bet)*sin(alp)-z*cos(gam)*sin(bet)*sin(alp)+x*sin(gam)*sin(bet)*sin(alp) z*sin(gam)*cos(alp) + x*cos(gam)*cos(alp) x"'=x* [sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)] + y* [cos(bet)*sin(alp) ] + z* [sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp) ] ---------------------------------------------------------- z"'=z" z"'=x* [- sin(gam)*cos(bet)] + y* [sin(bet)] + z* [cos(gam)*cos(bet)] ---------------------------------------------------------------- This really appear to be more complicated that what we had before. It appear to have much more multiplications than we had. But if we look on those resulting formulas we can see that all coefficients in square brackets can be calculated just once so that transformation of a point would look like: x"'=x*mx1 + y*my1 + z*mz1 y"'=x*mx2 + y*my2 + z*mz2 z"'=x*mx3 + y*my3 + z*mz3 It can be seen that it takes just 9 multiplication. Of course calculating all the coefficients takes 16 multiplications: mx1= sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp) my1= cos(bet)*sin(alp) mz1= sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp) mx2= sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp) my2= cos(bet)*cos(alp) mz2= -cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp) mx3= -sin(gam)*cos(bet) my3= sin(bet) mz3= cos(gam)*cos(bet) But if we have 100 points to turn. Original method would take 12*100=1200 multiplications. This one would take 9*100+16=916 since coefficients are calculated just once for all points to transform. One more consideration having to do with optimization, in most cases it makes sense to measure angles as integers. We don't usually care of fractions smaller then 1 degree. So we don't really need sin and cos for all the real number range. That's why we can calculate sin/cos just once for all 360 degrees before any rotations are performed and put the obtained values into arrays so that next time we would need any of them we wouldn't have to call a functions which probably uses a power series with few multiplications to calculate each. (Processors nowadays come with floating point units which can calculate sin/cos pretty fast but unlikely faster then single array access (this might become false though). By the way we don't really need 360 degrees. This number is not very convenient. We can go with full angle divided into 256 pseudo degrees. The reason we would need just one unsigned char to store the angle and whenever we go beyond 255 we simple wrap around to zero, saving this way a conditional statement we would need otherwise in the case of 360 degrees. #define T_RADS 40.7436642 /* pseudo grads into rads */ float T_mx1,T_mx2,T_mx3,T_my1,T_my2,T_my3,T_mz1,T_mz2,T_mz3; float T_sin[256],T_cos[256]; /* precalculated */ void T_init_math(void) { int i; for(i=0;i<256;i++) { T_sin[i]=sin(i/T_RADS); T_cos[i]=cos(i/T_RADS); } } Now finally to the code to perform 3D rotation. First function is building set of rotation coefficients. Note that T_mx1, T_mx2 etc. are global variables. The other function: T_rotation uses those coefficients, so they better be ready when later is called. void T_set_rotation_gam_bet_alp(int alp,int bet,int gam) { T_cosalp=T_cos[alp]; T_sinalp=T_sin[alp]; T_cosbet=T_cos[bet]; T_sinbet=T_sin[bet]; T_cosgam=T_cos[gam]; T_singam=T_sin[gam]; /* initializing */ T_mx1=T_singam*T_sinbet*T_sinalp + T_cosgam*T_cosalp; T_my1=T_cosbet*T_sinalp; T_mz1=T_singam*T_cosalp - T_cosgam*T_sinbet*T_sinalp; T_mx2=T_singam*T_sinbet*T_cosalp - T_cosgam*T_sinalp; T_my2=T_cosbet*T_cosalp; T_mz2=-T_cosgam*T_sinbet*T_cosalp - T_singam*T_sinalp; T_mx3=-T_singam*T_cosbet; T_my3=T_sinbet; T_mz3=T_cosgam*T_cosbet; /* calculating the matrix */ } void T_rotation(int *from,register int *to,int number) { register int i; register int xt,yt,zt; for(i=0;i| * | /| A (X,Z) | / | |/ | ^ * | | /|X' | / | | |/ | | *---+---+- - - - -> 0 | Z | | focus distance ->|---|<- pic 2.11 Let's consider 2-D case first: the viewer's eye is located at point 0 the distance between viewer's eye and the screen would be called - "focus". The problem is to determine at what point the ray going from point A into the viewer's eye would intersect the screen. We would have to plot the dot just there. This is very simple, just yet another case of solving similar triangles. Just from the fact that the angle at 0 is the same for both triangles, and if two angles are the same their tangents would be too. We have: X' X X*focus ----- = ----- => X'= ----------- focus Z Z Same considerations apply for Y dimension, together describing 3-D case: Y'=Y*focus/Z X'=X*focus/Z It is a bit of a mystery what is happening with Z coordinate. Basically after the perspective transformation we would be performing rendering onto the screen which is 2-D surface we would need X and Y but hardly Z, so we might as well have Z'=0. However when trying to render multiface objects we would need to know the depth (Z) of all polygons so that we can decide which are visible and which are obscured. So we might have Z'=Z that is just leave old value. Then again by doing transformation on X and Y and leaving Z unchanged we actually may allow for depth relation to change. That is a polygon obscuring another one in the original space, before projection, can actually come out behind or partly behind with this kind of transformation pic 2.21. ^ Z same | * * deapth | * / * / before | A / / // and | * / /* A after | /B / B | * * | before after pic 2.12 To avoid that we can apply some non linear transformation on Z, like Z'=C1+C2/Z, where C1, C2 are some constants. In the enclosed code only X and Y are being transformed. The implementation of all those things is that easy I won't even put a code example here. However, there are few possible questions: first: what should be the value for "focus distance" ? \ / \ / \ / \ / \ / \ / ---I\--------/I--- ---I\--------/I--- \ / \ / \/ \ / \ / \/ pic 2.13 pic 2.13 shows it's physical sense, focus basically determines the view angle. when it is smaller - angle is wider, when it is bigger angle is smaller. It is helpful to think of this distance as being measured in screen pixels, this way for 320 pixel display if focus is chosen to be 180 the view angle we are getting is 90'. Perspective transformation might produce at first, weird looking results, with distortions sometimes similar to that of wide-range photographic pictures. One has to experiment a bit to find what's looking best, values giving view angle somewhere between 75-85' would probably look not that terrible, but that depends of scene and screen geometries. clipping problem... ------------------- Where transforms point with Z equal to 0? since we are dividing by Z that would be infinity, or at least something producing divide error, in this computer-down-to-earth-world. The only way to deal with this is to guarantee that there would be no coordinates with such Zs. another thing, calculations for points with negative Z would produce negative coordinates, what we would see is objects (or even parts thereof) flipped over, but then again objects with negative coordinates are behind the viewing plane and effectively behind the viewer, so this is something we can't see and better make sure this is not being rendered. The way to do it is by applying 3-D clipping. The perspective transformation can also be represented as a matrix, however, let's think of all transformations discussed up to this moment a point would have to go through in order to be rendered. First we would rotate the world with respect to viewer orientation, then before we can apply perspective transformation we would have to do clipping and at least get rid of the vertices behind the viewing plane, and only then can we apply perspective transformation. Representing clipping in matrix form is not the most pleasant thing there is. And we've already discussed that matrices would be of use when we have several consecutive transformations so that all of them would be represented by one matrix. In this case however rotation and perspective are separated by the clipping. If we are sure that none of the coordinates would ever get behind the viewing plane, we don't actually need clipping, that's the instance to represent everything in the matrix form, otherwise it might not be a good idea neither in terms of code complicity nor performance. fixed point math. ----------------- In previous chapters in several places we didn't have another choice but to use floating point multiplications. (sin and cos are after all real-valued functions) However this is still fairly expensive operation, integer multiplications and divisions would usually be faster. However since we have to represent numbers with fractions it appeared we didn't have another choice but to use floating point operations. The question is, can we represent fractions using just integers? One might think there's a contradiction in the question itself, however we would see that, as in many other cases, we can use one mechanism to represent something else. This is the case here, We may use fixed point numbers representing them as integers, and with small amendments we can use regular integer additions multiplications etc, to perform fixed point operations. First of all how do we represent integers? +------+------+------+------+------+ | 1 | 0 | 0 | 0 | 1 | +------+------+------+------+------+ 4 3 2 1 0 2 =16 2 =8 2 =4 2 =2 2 =1 pic 2.14 With any counting system, digits in multi digit numbers would have different weight, so to speak. With decimal numbers this weight would be some power of ten, that is 102 is actually: 10**2*1 + 10**1*0 + 10**0*2 == 100+0+2 ==102 Same with binary numbers, the only difference, weight would be some power of 2. this way number on pic 2.13 would be: 2**4*1 + 2**3*0 + 2**2*0 + 2**1*0 + 2**0*1 == 16+0+0+0+1 == 17 Now, it is quite natural for us to place a decimal point into a decimal number, placing a binary point into a binary number should not seem more difficult. How do we interpret the decimal point? It is just that for the numbers to the right of the decimal point we pick negative power of ten. That is: 102.15 == 10**2*1 + 10**1*0 + 10**0*2 + 10**-1*1 + 10**-2*5 == 1 5 15 == 100 + 0 + 2 + ---- + ---- == 102 + ----- 10 100 100 In this representation the point separates negative powers from zero power and unlike numbers represented in exponential form - the floating point numbers, in this case the point never changes it's position. Nothing should stop us from extending the same scheme onto binary numbers: +------+------+------+------+------+ | 1 | 0 | 0 | 1 | 1 | +------+------+------+------+------+ 2 1 0 -1 -2 2 =4 2 =2 2 =1 2 =1/2 2 =1/4 pic 2.15 Using the very same technique we can see that number represented on the pic 2.15 can be considered also as: 2**2*1 + 2**1*0 + 2**0*0 + 2**-1*1 + 2**-2*1 == 4+0+0+1/2+1/4 == 4 + 3/4 Numbers in what range can we represent? Resorting once again to analogy with decimal numbers it can be seen that two digits to the right from the decimal point can cover the range of 0.99 in 0.01 steps. Numbers smaller then minimal 0.01 precision step just can't be represented without increasing number of digits after the decimal point. Very same thing is happening with binaries in the example at pic 2.15 since we have two binary digits- minimal number we can represent is 1/4 ( 0.01 (bin) == 1/4 (dec) ) and again numbers with higher precision can be represented only by increasing number of binary digits after the binary point. Now to the representation of negative numbers, there actually several ways to do that but effectively today one has prevailed, this is so called 2's compliment form. The idea is to negate the weight of the leftmost digit pic 2.16: +------+------+------+------+------+ | 1 | 1 | 1 | 1 | 1 | +------+------+------+------+------+ 4 3 2 1 0 -2 =-16 2 =8 2 =4 2 =2 2 =1 pic 2.16 the number represented above is considered to be: -2**4*1 + 2**3*1 + 2**2*1 + 2**1*1 + 2**0*1== -16+8+4+2+1 = -1 The advantage of this approach, we don't have to change addition and subtraction algorithm working for positive integers in order to accommodate 2's compliment numbers. Why is that? just by the virtue of natural overwrap of integer numbers the way they are represented in the computers. If we add 1 to the maximum number we can represent we would obtain a carry from the most significant (leftmost) bit, which is ignored, and the value of exactly 0. -1 in signed representation is exactly the maximum possible value to represent in unsigned representation, -1+1==0 indeed. (One can check this to be true for all numbers and by more formal means, but since this is not the point I would ignore it) Again, addition and subtraction don't have to be changed to accommodate 2's compliment negative numbers (multiplication however should, that's why there are instructions for both signed and unsigned multiplications in most computers) Since the leftmost digit in 2's compliment representation carries negative weight and that's exactly the one with highest power the minimum negative number possible to represent in the example above would be 10000 (bin) = -16 all the other digits have positive weights so maximum possible positive number would be 01111 (bin) = 15 but this slight asymmetry is not a problem in the majority of cases. However, values of sin and cosin functions are between 1 and -1 so to represent them we might choose format with just one integer field: +------+------+------+------+------+ | 0 | 1 | 1 | 1 | 1 | +------+------+------+------+------+ 0 -1 -2 -3 -4 -2 =-1 2 =1/2 2 =1/4 2 =1/8 2 =1/16 pic 2.17 but, because of the asymmetry described above there would be a value for -1 (1000) but there would be no pure 1, just it's approximation (01111): (pic 2.17) 1/2+1/4+1/8+1/16, Then again, for most graphics applications when we have, say, 15 fields representing fractions this would be close enough and won't cause a lot of problems. As you can see, the very same rules work for 2's compliment whether with fixed point or without. multiplication of fixed point numbers. -------------------------------------- Let's consider what is happening when we are multiplying two decimal numbers, for example we have to multiply an integer by a number with a decimal point and we would need just an integer as a result: 1.52 * 11 ----- 1 52 + 15 2 ----- 16.72 Actual result of multiplications has as many digits after the decimal point as there were in both arguments. Since we need just an integer as a result we would discard numbers after the point effectively shift digits to the right by 2 in this case. We can do the same with fixed point binary numbers, it means that we have to readjust precision after the multiplication. Say, if we are multiplying an integer by a number having 8bit considered to be after the binary point the result would also 8bit after the point. If it is just the integer part we are interested in, we have to shift the result right 8 bit destroying all fractional information. Similar things are happening with division except that if we want to divide two integers and obtain a fixed point result the argument to be divided has to be shifted left- added fixed point fields, so that effectively we are dividing a fixed point number by an integer. How to implement all that? Well, first of all we have to decide what kind of values and what operations on them we would need. Sometimes it is nice just to have general purpose fixed point math library, on the other hand we may choose several different formats of fixed point numbers one for every special place. Another thing we have to consider is: would we be able to represent all the operations using just high level C constructs or it would be necessary to descend into assembly instructions. The reason we might consider later is: in most assemblies the result of multiplication has twice as many bits as each of arguments, and since we would occupy some bits with fractions we really need all bits we can get. Another thing, the result of multiplication would often be located in two registers, So we can organize operations the way so that instead of adjusting result by shifting we would just take it from the register carrying higher bits, -effectively doing zero-cost shifting right by the bit length of the register. On the other hand, do we really want to compromise code portability by coding in assembly? Obviously, there are valid reasons to answer that both ways. In particular case of this text and associated code we are interested in good portability so let's try to stay within straight C boundaries. Where would we be using fixed point? first of all - in 3-D transformations, to keep rotation matrix coefficients. We would need two kinds of multiplication: In calculation of the coefficients: fixed * fixed producing same precision fixed and in actual coordinate transformation fixed * int producing int. lets define fixed to be: typedef signed_32_bit fixed; /* better be 32bit machine */ #define T_P 10 /* fixed math precision */ the T_P would be the number of bits we are allocating to keep fractions. Let's make a little function to transform the floating point values in the range [-1,1] into fixed point values: fixed TI_float2fixed(float a) { fixed res=0; int i,dv,sign; if((sign=a<0.0)==1) a=-a; for(i=1,dv=2;i<=T_P;i++,dv*=2) if(a>=1.0/dv) { a-=1.0/dv; res|=(0x1<<(T_P-i)); }; if(sign==1) res=-res; return(res); } In the function above we are iteratively subtracting powers of two from the argument and using the result to set the bits in the fixed point number. We would use this function to precalculate the array of sin/cos values. When calculating the rotational coefficients we would multiply one fixed by another and obtain same precision fixed result. again, actual result would have twice as many fractional bits as each of the argument so same precision result would be (a*b)>>T_P where 'a' and 'b' fixed point values and T_P number of bits after the binary point in each. This way code would look something like: ... ... T_wx1=((singam*((sinbet*sinalp)>>T_P))>>T_P) + ((cosgam*cosalp)>>T_P); T_wy1=((cosbet*sinalp)>>T_P); T_wz1=((singam*cosalp)>>T_P) - ((cosgam*((sinbet*sinalp)>>T_P))>>T_P); ... ... In rotation function ints would have to be multiplied by a fixed coefficient, the result would have T_P fractional bits, since we are interested in just integer part same rule as above applies: ... ... *to++=((T_wx1*xt)>>T_P) + ((T_wy1*yt)>>T_P) + ((T_wz1*zt)>>T_P); *to++=((T_wx2*xt)>>T_P) + ((T_wy2*yt)>>T_P) + ((T_wz2*zt)>>T_P); *to++=((T_wx3*xt)>>T_P) + ((T_wy3*yt)>>T_P) + ((T_wz3*zt)>>T_P); ... ... That's about it. The only thing, we really have to be careful with the range of numbers the operations would work for. If we are working with 32bit numbers and would choose 16 of them for fractional part and one for sign bit then only 15bit would carry the integer part. after multiplication of two fixed numbers since result has as many fractional bits as in both arguments all 32 bits would actually carry fractions. The integer part would be gone in the dark realm of left-carry, and there's no way in standard straight C to get it out of there. (when using assembly we don't have this problem since result of multiplication has twice bits of each argument, so we really can get trough to the integer part, however, most C compilers define the result of multiplication to be of the same bit size as arguments int*int==int, some compilers try to be smarter then that but we really have to consider worst-case scenario) And finally is it really true that integer multiplication plus shifts is faster then floating multiplication? What follows are results of some very approximate tests: sparc: floating faster fixed about 1.3 times 68040: floating faster fixed about 1.5 times rs4000: floating faster fixed about 1.1 times rs6000: fixed faster floating about 1.1 times i386sx: fixed faster floating about 5.1 times floating in the test was: float a,b; a*b; fixed in the test was: int a,b; (a*b)>>15; As one can see processors with fast floating point units nowadays do floating operations faster then integer once. Then again in terms of portability integer multiplication would always be fast enough whereas floating multiplication especially on processors without FPUs (i386sx) might get really slow. The decision whether to implement fixed point, this way, might get a bit tough to make. In the provided library's source transformations are implemented both ways but with the same function prototypes, so that one can decide which implementation is more suitable to link in. * * *